A. Polynomial-based interpolation

Image interpolation aims to produce a high-resolution image by upsampling the low-resolution image. As explained in the introduction section, interpolation algorithms often assume that the observed LR image is a direct downsampled version of the HR image. Hence, the de-aliasing ability during the upsampling process is important, i.e. the recovery of the high-frequency signal from the aliased low-frequency signal [1].

For real-time applications, conventional polynomial-based interpolation methods such as bilinear and bicubic interpolation are often used due to their computational simplicity [2]–​[5]. For example, the bicubic convolution interpolator requires only several arithmetic operations per pixel, such that real-time processing can easily be achieved [4]. The basic idea is to model the image signal by a low-order polynomial function (using some observed signals). However, polynomial functions are not good at modeling the signal's discontinuities (e.g. edges). Hence, the conventional polynomial-based interpolation methods often produce annoying artifacts such as aliasing, blur, halo, etc. around the edges. To resolve this problem, some adaptive polynomial-based and step function-based interpolation methods were proposed [6]–​[9].

B. Edge-directed interpolation

B3. New edge-directed interpolation (NEDI)

New edge-directed interpolation (NEDI) uses the FIR Wiener filter, equivalently, the linear minimum mean squares error estimator [19] for linear prediction. Let us briefly describe the formulation of NEDI as an introduction. The 4-order linear estimation model is given by

Yj=∑k=14Ak⋅YiΔk+εi,fori=1,2,…,P(2.1)

View SourceY_{j}=\sum_{k=1}^{{4}}A_{k}\cdot Y_{i\Delta k}+\varepsilon_{i},\quad for\quad i=1,2,\ldots,P\eqno\hbox{(2.1)}where εi is the estimation error, P is the number of data point samples, Ak is the model parameters and each available data point sample (Yi) has four neighboring data points, YiΔk. Figure 1a shows the spatial positions of YiΔk and Yi. The model parameters Ak can be found by using the least squares estimation:

{A^k}=argmin{A^k}∑i=1P[Yi−∑k=14Ak⋅YiΔk]2(2.2)

View Source\{\hat{A}_{k}\}=\mathop{\arg\min}_{\{\hat{A}_{k}\}}\sum_{i=1}^{P}\left[Y_{i}-\sum_{k=1}^{4}A_{k}\cdot Y_{i\Delta k}\right]^{2} \eqno\hbox{(2.2)}where the matrix form of (2.2) is given by

A^=argminA∥Y−YAA∥22(2.3)

View Source\hat{{\bf A}}=\mathop{\arg\min}_{\bf A}\Vert {\bf Y}-{\bf Y}_{{\bf A}}{\bf A}\Vert_{2}^{2}\eqno\hbox{(2.3)}and the matrices are defined as

Y={Yi}T,YA={YiΔk}T,A={Ak}T(2.4)

View Source{\bf Y}=\{Y_{i}\}^{T},{\bf Y}_{{\bf A}}=\{Y_{i\Delta k}\}^{T},{\bf A}=\{A_{k}\}^{T} \eqno\hbox{(2.4)}

The sizes of matrices Y,YA and A are P×1,P×4 and 4xl respectively. The close form solution of (2.3) is given as

A^=(YTAYA)−1YTAY(2.5)

View Source\hat{{\bf A}}=({\bf Y}_{{\bf A}}^{T}{\bf Y}_{{\bf A}})^{-1}{\bf Y}_{{\bf A}}^{T}{\bf Y} \eqno\hbox{(2.5)}where A^ is called the ordinary least squares (OLS) estimator. Due to “geometry duality” [20], the missing data point X can be interpolated by its four neighboring data points, Xk, as follows:

X=∑k=14A^k⋅Xk(2.6)

View SourceX=\sum_{k=1}^{4}\hat{A}_{k}\cdot X_{k} \eqno\hbox{(2.6)}

For the missing data point between two available LR data points (e.g. the missing data point between X1 and X2 in Figure 1b), its value is obtained by rotating the spatial positions of the neighbors and missing data points by 45 degree with a scaling factor of 1/2–√. More details can be found in [19]. Note that the linear interpolation method described in this paper can be applied to high activity areas only, which can be identified by using a local variance larger than 8, as used in NEDI. The number of samples is P=64, i.e. using an 8×8 window.

Fig. 1 Graphical illustration of spatial positions of LR and HR pixels.

Show All

NEDI uses the LR image to estimate the HR covariances for the FIR Wiener filter, as shown in (2.5) and (2.6). There are many algorithms which are based on the idea of NEDI [21]–​[30]. In [21], a six-order filter design was proposed to avoid error accumulation in the original two steps design. In [22], an eight-order filter was also proposed to include the vertical-horizontal correlation. The filter weights of the FIR Wiener filter were regularized for stability [23]. In [31], instead of using the LR image to estimate the HR covariances; the relevant HR image pairs were searched from an offline dictionary to estimate the HR covariances. Let us clarify this filtering process with the following over-simplified example.

Example:

A data sequence is given by {Y(0),Y(1),Y(2),Y(3), Y(4),Y(5),Y(6),Y(7),Y(8),Y(9)}. We have to find Y(10) using the FIR Wiener filter. Let also the length of the filter be 2, i.e. N=2 Hence the observed sample data point Y and their corresponding interpolating points, YA, are as shown below

Y=[Y(2)Y(3)Y(4)Y(5)Y(6)Y(7)Y(8)Y(9)]TYA=[Y(0)Y(1)Y(1)Y(2)Y(2)Y(3)Y(3)Y(4)Y(4)Y(5)Y(5)Y(6)Y(6)Y(7)Y(7)Y(8)]T

View Source\eqalignno{ & {\bf Y}=[\matrix{Y(2)& Y(3) & Y(4) & Y(5) & Y(6) & Y(7) & Y(8) & Y(9)]^{{\bf T}}}\cr & {\bf Y}_{{\bf A}}=\left[\matrix{Y(0) & Y(1) & Y(2) & Y(3) & Y(4) & Y(5) & Y(6) & Y(7)\cr Y(1) & Y(2) & Y(3) & Y(4) & Y(5) & Y(6) & Y(7) & Y(8)}\right]^{T}}This means that {Y(9),Y(8),…Y(2)} are used as sample points to look for the required statistics. The interpolation points for Y(9) are Y(7) and Y(8); for Y(8) are Y(6) and Y(7); … for Y(2) are Ye(0) and Y(1). Let us assume further that the cross-correlation matrix YATY=[18 11]T and the auto-correlation matrix YTAYA=[13885]. The filter weights can be estimated by equation (2.5), as follows: A=[A1A2]=(YTAYA)−1YTAY=[13885]−1[1811]=[2−1]. Hence, the data point Y(10) can be estimated by (eqn.2.6), i.e.

Y(10)=A1Y(8)+A2Y(9)=2Y(8)−Y(9).

View SourceY(10)=A_{1}Y(8)+A_{2}Y(9)=2Y(8)-Y(9).

B4. Soft-decision interpolation (SAI)

The soft-decision adaptive interpolation (SAI) [24] was proposed to interpolate a block of pixels at one time using the idea of NEDI. It has the benefit of using the block-based estimation, by constraining the statistical consistency within the block region, which comprises of both observed LR and unobserved HR pixels. Within a local window as shown in Figure 2, the observed LR pixels Y in the LR image are used to estimate the HR pixels X in the HR image. The sizes of X and Y in the SAI are 12×1 and 21×1 respectively. Let us formulate the SAI as a maximum a posterior (MAP) estimation problem, as follows

argmaxx p(X/Y)=argmaxxlog[p(Y/X)p(X)]=argminx(∥X−YxA∥22+∥Y~−XYA∥22+λ∥X~−XxB∥22)(2.7)

View Source\eqalignno{ & \mathop{\arg\max}_{\rm x}\ p({\bf X}/{\bf Y})=\mathop{\arg\max}_{\rm x} \log [p({\bf Y}/{\bf X})p({\bf X})]\cr & =\mathop{\arg\min}_{\rm x}\left(\Vert {\bf X}-{\bf Y}_{{\bf x}}{\bf A}\Vert_{2}^{2}+\Vert\tilde{{\bf Y}}-{\bf X}_{{\bf Y}}{\bf A}\Vert_{2}^{2}+\lambda\Vert\tilde{{\bf X}}-{\bf X}_{{\bf x}}{\bf B}\Vert_{2}^{2}\right) &\hbox{(2.7)}}where the likelihood and prior are assumed to follow the Gaussian distributions, elements of A and B are the model parameters (which are estimated using the NEDI method), YX are diagonal neighbors of X, Y; are the centermost five elements of Y (bounded by dotted cross in Figure 2), XY are diagonal neighbors of Y~,X~ is a 4×1 vector representing the centermost four elements of X (bounded by dotted square in Figure 2), and XX are horizontal-vertical neighbors of X~. Hence, all the elements of X and Y are defined in the posterior. λ is the regularization factor which is recommended to be 0.5 [24]. In (2.7), the SAI uses three auto-regressive models to model the image signals. Such a formulation constrains the statistical consistency within a local region, which is usually locally stationary. Hence, this method is called the soft-decision adaptive interpolation. To solve (2.7), we can use a more compact form of the argument, as follows

argminx(CX− DY)T(CX−DY)(2.8)

View Source\mathop{\arg\min}_{\rm x} ({\bf CX}-\ {\bf DY})^{T} ({\bf CX}-{\bf DY})\eqno\hbox{(2.8)}which remaps the elements into a very compact form. Matrices C and D are defined as, C=[I12C1λC2]T and D=[D1D204×21]T (see reference [24] for details). The argument in (2.8) can be solved by differentiating the cost function with respect to X to obtain the following closed-form solution

X=(CTC)−1CTDY(2.9)

View Source{\bf X}=({\bf C}^{T}{\bf C})^{-1}{\bf C}^{T} {\bf DY}\eqno\hbox{(2.9)}

Fig. 2 Graphical illustration of 21 LR and 12 HR pixels within the local window used by SAI. Left figure shows one element of Xi∈X and Yi∈Y, and their diagonal neighbors

Show All

B5. Robust soft-decision interpolation using weighted least squares (RSAI)

It is well known that the least squares estimation is not robust to outliers, hence the weighted least squares was proposed to improve the accuracy and robustness of the SAl [25]–​[26]. The robust SAI (RSAI) incorporates the weights to all residuals in the cost function of the SAI, as follows

argminx(∥W1(X−YxA)∥22+∥W2(Y~−XYA)∥22+λ∥W3(X~−XxB)∥22)(2.10)

View Source\mathop{\arg\min}_{\rm x}\left(\matrix{\Vert {\bf W}_{1}({\bf X}-{\bf Y}_{{\rm x}}{\bf A})\Vert_{2}^{2}+\Vert {\bf W}_{2}(\tilde{{\bf Y}}-{\bf X}_{{\bf Y}}{\bf A})\Vert_{2}^{2}\cr \hfill +\lambda\Vert {\bf W}_{3}(\tilde{{\bf X}}-{\bf X}_{{\bf x}}{\bf B})\Vert_{2}^{2}}\right)\eqno\hbox{(2.10)}where Wi contains the weighting parameters for residuals in (X−YXA),(Y~−XYA) and (X~−XXB) respectively. To solve (2.10), equations similar to (2.8)–​(2.9) can be used to derive the cost function and the resultant closed-form solution. As verified in [25], there exists severe outliers due to the mismatch of “geometric duality”, which is a fundamental assumption of estimating the HR parameters using the LR image. Hence, using the weighted least squares to adaptively weight the residuals for both HR parameters and HR pixels estimation results in a significant improvement of PSNR and SSIM, as shown in table 1 and table 2. As a matter of fact that the RSAI is currently one of the best performers for image interpolation.

Table I PSNR (dB) of the different image interpolation algorithms

Table II SSIM [79] of the different image interpolation algorithms

B6. Bilateral soft-decision interpolation for real-time applications (BSAI)

The SAI and RSAI are able to produce the best quality, but the computational cost is high. To largely reduce the computation, the least squares estimation should be avoided or reduced to single parameter estimation. Hence, the bilateral soft-decision interpolation (BSAI) was proposed to use the bilateral filter for replacing the least squares estimation for the model parameters and to reduce the weighted least squares soft-decision estimation to a single parameter estimation [17]. The BSAI uses bilateral filter weights Ak to replace the least squares parameters and adopts the following cost function (modified from SAI and RSAI) for estimating a HR pixel at one time (instead of a block of HR pixels at one time in the SAI and RSAI)

argminX⎧⎩⎨(X−∑k=03AkXk)2+∑k=03Uk(Xk−∑i=03A,Xki)2⎫⎭⎬(2.11)

View Source\mathop{\arg\min}_{X}\left\{\left(X-\sum_{k=0}^{3}A_{k}X_{k}\right)^{2}+\sum_{k=0}^{3}U_{k}\left(X_{k}-\sum_{i=0}^{3}A,X_{ki}\right)^{2}\right\} \eqno\hbox{(2.11)}where Xk 's are the neighbors of X which is the HR pixel to be interpolated), Xki 's are the neighbors of Xk. Figure 3 shows the spatial configurations of Xk and X, and X0i and X0 as examples. In (2.11), the first term constrains the result to approach the interpolated result using the bilateral filter, and the second term constrains the interpolated result to be continuous with its neighbors (soft-decision estimation). The weight Uk is defined by the bilateral filter weights (since the continuity property of edge depends on the edge orientation, which is exploited by the bilateral filter weight) added with a constant for stabilization. To solve (2.11), we differentiate the cost function with respect to the variable X and obtain the following closed-form solution,

X=∑k=03AkXk+∑k=0;l=3−k3UkAl(Xk−∑k=0;l≠3−k3AiXki)1+U0A23+U1A22+U2A21+U3A20(2.12)

View SourceX={\sum\limits_{k=0}^{3}A_{k}X_{k}+\sum\limits_{k=0;l=3-k}^{3}U_{k}A_{l}\left(X_{k}-\sum\limits_{{k=0;l\neq 3-k}}^{3}A_{i}X_{ki}\right)\over 1+U_{0}A_{3}^{2}+U_{1}A_{2}^{2}+U_{2}A_{1}^{2}+U_{3}A_{0}^{2}} \eqno\hbox{(2.12)}which is a very compact and efficient equation. In real-time implementation using C++ codes, BSAI requires 0.062 second to interpolate a 384×256 image to double its size. The major advantage of BSAI is that it is the near-best performer for the image interpolation (as shown in table 1 to table 3. but its computational cost is sufficiently low for real-time applications.

Fig. 3 Spatial configurations of XkX (left), and X0iX0 (right).

Show All

B7. Recent trends af image interpolation applications

Recently, image interpolation does benefit from the development of sparse representation. A family of linear estimators corresponding to different priors are mixed using the sparse representation to give the final estimates for image interpolation [34]. Its performance is on a par with the SAI [24]. Obviously, the sparse representation, which has shown to be successful in compressive sensing, denoising, restoration and super-resolution, etc, is a fruitful research direction of image interpolation.

Another trend of the image interpolation application is the real-time upsampling of a LR video sequence for future very high-definition TVs, such as videos with 4K resolution. The HR video sequence can be obtained by fusing several lowresolution (LR) frames into one high-resolution (HR) frame. For such real-time applications, nonlocal means (weighted sum filter) can be directly applied due to its simplicity in computing the filter weights and its high performance. The nonlocal means can be used with a linear motion model to better estimate the filter weights, which is optimized for the upsampling [35]. After fast upsampling, some simple restoration techniques can be applied to deblur the video. In the future, there should be more of these fast algorithms which can perform upsampling in real-time, of which the realtime requirement is the major difficulty for their development.

Recently view synthesis attracts much attention in the image processing community. In 3D videos and multi-view synthesis, the major issue is to fill the newly exposed area (hole region). The hole filling in view synthesis resembles the video (multi-frame) interpolation scenario. Hence, the multiframe interpolation techniques such as nonlocal means can be adopted to fill the hole. Specifically, the nonlocal means can be modified to cope with irregular hole sizes and the extra depth information [36]. Similar to eqn. (2.6), the nonlocal means for hole filling can be defined as

x=∑i∈Nwixi(2.13)

View Sourcex=\sum_{i\in N}w_{i}x_{i}\eqno\hbox{(2.13)}where x is a pixel's intensity inside the hole region, xi are neighbors of x,N is the set of neighbors and wi is the weight of the neighbor xi. A new closeness term (which comprises of the intensity and depth difference) was proposed to determine the weights of non local means [36], as follows

c(x,xi)=∑j∈Wx;k∈Wx(i)[x(j)−xi(k)]2[dx(j)−dx(i)(k)]2(2.14)

View Sourcec(x,x_{i})=\sum_{j\in W_{x};k\in W_{x(i)}}[x(j)-x_{i}(k)]^{2}[d_{x}(j)-d_{x(i)}(k)]^{2}\eqno\hbox{(2.14)}where x(j) and dx(j) represent the nearby intensity and depth values of x within the local window Wx, while xi(k) and dx(i)(k) represent the nearby intensity and depth values of xi within the local window Wx(i). After some further development, the weight of the modified nonlocal means for hole filling is defined as

wi=exp(−c(x,xi)/σc)⋅exp(−d(x,xi)/σd)∑i∈Nexp(−c(x,xi)/σc)⋅exp(−d(x,xi)/σd)(2.15)

View Sourcew_{i}={\exp(-c(x,x_{i})/\sigma_{c})\cdot\exp(-d(x,x_{i})/\sigma_{d})\over \sum\limits_{i\in N}\exp(-c(x,x_{i})/\sigma_{c})\cdot\exp(-d(x,x_{i})/\sigma_{d})} \eqno\hbox{(2.15)}where d(.,.) measures the geometric distance between two pixels and σ stands for the variance. Experimental results show that this new approach outperforms the conventional spatial and temporal approaches for hole filling during the view synthesis. However, we believe that more interpolation techniques can be novelty designed suitable for view synthesis, which is a good future direction for image interpolation.

C. PSNR and SSIM comparisons af several image interpolation algorithms

Let us compare the subjective and objective quality of some image interpolation algorithms as described in this paper. Table I and Table II show the PSNR and SSIM values of various image interpolation algorithms [4], [15], [17], [19], [24], [25] using 24 natural images from Kodak. We have also considered the execution time of these image interpolation algorithms using C++ codes as shown in Table III (see also figure 5). It is shown that the RSAI [25] achieves the highest average PSNR and SSIM but its execution time is longer than the second and the third best performers, i.e. SAI [24] and BSAI [17]. BSAI requires much less computation compared to that of the SAI but its quality is closed to the SAI For realtime applications, BSAI is the best choice; and for offline applications, RSAI is the best choice.

Table III Execution time (second) of image interpolation algorithms using C++ codes (* is the estimated time)

D. Classification ofimage interpolation algorithms

In this section, we are going to classify broadly image interpolation algorithms, as shown in figure 4. They are the polynomial-based approaches which are fast and can be adaptive to local statistics, and the edge-directed approaches which directly address the edge reconstruction criterion. Edge-directed methods can intuitively interpolate along one edge orientation, or fuse several estimates of different edge orientations, or minimize the linear mean squares error, for example, using the FIR Wiener filter, whereas the soft-decision interpolation can also be applied.

Fig. 4 Classifications of the image interpolation algorithms

Show All

Fig. 5 Plot of PSNR and SSIM [79] against the execution time of image interpolation algorithms

Show All

A. Single-image and multi-frame super-resolution

Super-resolution (SR) aims to produce a high-resolution (HR) image using one or several observed low-resolution (LR) images by upsampling, deblurring and denoising. For multi-frame SR, there can involve registration and fusion processes. It is interesting to point out that some SR algorithms do not involve the denoising process, or some interpolation algorithms are also referred to as superresolution algorithms. Generally speaking, the superresolution methods can be classified into single-image superresolution (only one LR image is observed) [37]–​[54] and multi-frame super-resolution [55]–​[71]. The multi-frame superresolution for a video sequence can moreover use a recursive estimation approach of video frames, which is called the dynamic super-resolution [55]–​[58]. Dynamic super-resolution makes use of the previously reconstructed HR frame to estimate the current HR frame.

B. Reconstruction-based and learning-based superresolution

There are two major categories of SR algorithms: they are the reconstruction-based and learning-based algorithms. Reconstruction-based algorithms [37]–​[39], [55]–​[65], [67]–​[71] use the data fidelity function, usually with the prior knowledge to regularize the ill-posed solution. Gradients (edges) are the main features to be used as the prior knowledge. Some more general prior knowledge using nonlocal self-similarity was also proposed [62]. Learning-based algorithms [40]–​[54], [66] moreover utilize the information from training images. Generally, when a LR image is observed, a patch within that LR image is extracted and searched within the trained dictionary for estimating a suitable HR patch that reconstructs the HR image. Some recent investigations [41]–​[43] make use of the sparse representation for estimating the HR pair and training the dictionary. By combining the use of prior knowledge and dictionary training, some unified frameworks [42]–​[44] for single-image SR were proposed to further regularize the estimated solution from the training sets.

C. Generic and face super-resolution

For the learning-based algorithms, they are most suitable for specific applications, such as face super-resolution. Since the face SR is crucial in applications like face recognition, surveillance, etc, a number of super-resolution methods were proposed for hallucinating the faces [48]–​[54]. These are training-based methods which make use of the common characteristics of human faces (e.g. eigenface) to design the formulation for dictionary training and HR image reconstruction.

D. Blind and non-blind super-resolution

The super-resolution methods can also be classified as blind and non-blind methods. The blind methods treat the point spread function (PSF) that represents the blur, and the registration parameters as variables to be estimated simultaneously with the HR image [67]–​[71]. The Bayesian or MAP framework is widely used to jointly estimate the variables. Non-blind methods are still widely used because the PSF can be approximated, can be set according to the user preference, or is known due to the knowledge of the camera. The registration parameters can also be separately estimated. Usually, the non-blind methods apply the regularizations to make them robust to inaccurate estimations of PSF, registration parameters, etc. The non-blind methods are widely used due to its simplicity in formulations, which are also easier to be parallelized for the practical applications.

E. FIR Wiener filter for super-resolution

The finite-impulse response (FIR) Wiener filter, or equivalently, the linear minimum mean squares error (LMMSE) estimator, can be applied to perform superresolution reconstruction. The block-based FIR Wiener filters were proposed for multi-frame SR [65]–​[66]. In [65], a wide-sense stationary correlation function based on the geometric distance between pixels was proposed to estimate the covariances for the FIR Wiener filter, and elegant results were reported. The partition filters partitions an image into blocks for applying the FIR Wiener filter, where the filter weights are learned from an offline dictionary, which is retrieved during the online estimation [66].

Similarly, the Gaussian process regression (GPR) can be applied for super-resolution reconstruction. GPR provides a sophisticated Bayesian framework to estimate the HR image and the hyper-parameters of the process (e.g. noise variance, and variance of the correlation function) using a pilot HR image, which can be obtained by the bicubic interpolation [37]. The nonlocal means has been proposed as the correlation function in GPR [37]. Using the nonlocal means as the correlation function, the iterative scheme of FIR Wiener filter can alternatively update the estimated covariances and HR image, which can address the disadvantage of inaccurate pilot HR image [47].

Let us briefly introduce the iterative Wiener filter (IWF) algorithm [47] which currently produces the best PSNR and SSIM results among some state-of-the-art algorithms, as shown in Table IV and Table V. The formulation of the FIR Wiener filter which minimizes the linear mean squares error is given by

View Source{\bf W}_{i}={\bf R}_{i}^{-1}{\bf p}_{i} \eqno\hbox{(3.1)}

y^i=WTixi(3.2)

View Source\hat{{\bf y}}_{i}={\bf W}_{i}^{T}{\bf x}_{i} \eqno\hbox{(3.2)}

Wi=[cmaxE{XiXTi/cmax}]−1cmaxE{xiyTi/cmax}=[E{xixTi/cmax}]−1E{xiyTi/cmax}(3.3)

View Source\eqalignno{ {\bf W}_{i} & =\left[c_{\max}E\{{\bf X}_{i}{\bf X}_{i}^{T}/c_{\max}\}\right]^{-1}c_{\max}E\{{\bf x}_{i}{\bf y}_{i}^{T}/c_{\max}\}\cr & =\left[E\{{\bf x}_{i}{\bf x}_{i}^{T}/c_{\max}\}\right]^{-1}E\{{\bf x}_{i}{\bf y}_{i}^{T}/c_{\max}\}&\hbox{(3.3)}}

E{pjpk/cmax}≈exp(−E(pj−pk)2/σ)(3.4)

View SourceE\{p_{j}p_{k}/c_{\max}\}\approx\exp(-E(p_{j}-p_{k})^{2}/\sigma) \eqno\hbox{(3.4)}

Algorithm 1 The iterative scheme of the FIR Wiener filter for SRR

Fig. 6 An example illustration of the observation vector xi and the unobserved vector yi when the magnification factor q=3 and number of elements in xi is 16.

Show All

Table IV PSNR (dB) of the estimated images using different algorithms

Table V SSIM [79] of the estimated images using different algorithms

Example:

A data sequence is given by {x(0),x(1),x(3),x(4)}, where x(2) is a missing data point between x(1) and x(3). We have to find x(2) using the Iterative Wiener filter. Let us define the observation vector xi and unobserved vector yi as below

xj=[x(1)x(3)]T and yi=x(2)

View Source{\bf x}_{j}=[x(1)\quad x(3)]^{T}\ {\rm and}\ {\bf y}_{i}=x(2)

Let us initialize the unobserved vector as the average of the nearest two data points as follows

View Source{\bf y}_{i}^{(0)}=(x(1)+x(3))/2

Let us compute the auto-correlation matrix and cross-correlation matrix using the correlation function in (3.4) with the nearest three data points for the expectation E(.):

First iteration: \

View Source\eqalignno{ {\bf R}_{i}^{(0)} & =\left[\matrix{ {\rm exp}(-E(x(1)-x(1))^{2}/\sigma) & {\rm exp}(-E(x(1)-x(3))^{2}/\sigma)\cr {\rm exp}(-E(x(3)-x(1))^{2}/\sigma) & {\rm exp}(-E(x(3)-x(3))^{2}/\sigma)}\right]\cr & =\left[\matrix{ {\rm exp}(-(0)/\sigma) & {\rm exp}(-E(x(1)-x(3))^{2}/\sigma)\cr {\rm exp}(-E(x(1)-x(3))^{2}/\sigma) & {\rm exp}(-(0)/\sigma)}\right]}

exp(−E(x(1)−x(3))2/σ)=exp(−E(x(3)−x(1))2/σ)=exp(−[(x(0)−y(0)j)2+(x(1)−x(3))2+(y(0)i−x(4))2/σ)

View Source\eqalignno{ & \exp(-E(x(1)-x(3))^{2}/\sigma)=\exp(-E(x(3)-x(1))^{2}/\sigma)\cr & =\exp(-[(x(0)-{\bf y}_{j}^{(0)})^{2}+(x(1)-x(3))^{2}+({\bf y}_{i}^{(0)}-x(4))^{2}/\sigma)}

R(0)i=[10.50.51]

View Source{\bf R}_{i}^{(0)}=\left[\matrix{1 & 0.5\cr 0.5 & 1}\right]

The cross correlation matrix is given by

View Source\eqalignno{ {\bf P}_{i}^{(0)} & =\left[\matrix{ {\rm exp}(-E(x(1)-x(2))^{2}/\sigma)\cr {\rm exp}(-E(x(3)-x(2))^{2}/\sigma)}\right]\cr & =\left[\matrix{ {\rm exp}(-[(x(0)-x(1))^{2}+(x(1)-{\bf y}_{j}^{(0)})^{2}+({\bf y}_{j}^{(0)}-x(3))^{2}/\sigma)\cr {\rm exp}(-[({\bf y}_{j}^{(0)}-x({\bf 1}))^{2}+(x(3)-{\bf y}_{i}^{(0)})^{2}+(x(4)-x(3))^{2}/\sigma)}\right]}